Mackey-Glass equation are a set of delayed differential equations describing the temporal behaviour of different physiological signal, for example, the relative quantity of mature blood cells over time.
The equations are defined as:
$$ \frac{dP(t)}{dt} = \frac{a P(t - \tau)}{1 + P(t - \tau)^n} - bP(t) $$where $a = 0.2$, $b = 0.1$, $n = 10$, and the time delay $\tau = 17$. $\tau$ controls the chaotic behaviour of the equations (the higher it is, the more chaotic the timeserie becomes. $\tau=17$ already gives good chaotic results.)
from reservoirpy.datasets.chaos import mackey_glass
timesteps = 25000
tau = 17
X, t = mackey_glass(timesteps, tau=tau)
# rescale between -1 and 1
X = 2 * (X - X.min()) / (X.max() - X.min()) - 1
plot_mackey_glass(X, t, 500, tau)
forecast = 10 # for now, predict 10 steps ahead
(X_train, y_train), (X_test, y_test) = split_timeserie_for_task1(forecast)
sample = 500
fig = plt.figure(figsize=(15, 5))
plt.plot(X_train[:sample], label="Training data")
plt.plot(y_train[:sample], label="True prediction")
plt.legend()
plt.show()
units = 100
leak_rate = 0.3
spectral_radius = 1.25
input_scaling = 1.0
density = 0.1
input_connectivity = 0.2
regularization = 1e-8
seed = 1234
im = plt.imread("./static/esn.png")
plt.figure(figsize=(15, 15)); plt.imshow(im); plt.axis('off'); plt.show()
Win = mat_gen.generate_input_weights(units, 1, input_scaling=input_scaling,
proba=input_connectivity, input_bias=True,
seed=seed)
W = mat_gen.generate_internal_weights(units, spectral_radius=spectral_radius,
proba=density, seed=seed)
reservoir = ESN(leak_rate, W, Win, ridge=regularization)
states = reservoir.train([X_train.reshape(-1, 1)], [y_train.reshape(-1, 1)], verbose=True)
y_pred, states1 = reservoir.run([X_test.reshape(-1, 1)], init_state=states[0][-1], verbose=True)
y_pred = y_pred[0].flatten()
states1 = states1[0]
Computing states: 100%|██████████| 1/1 [00:00<00:00, 2046.00it/s]
Training on 1 inputs (20000 steps)-- wash: 0 steps
Computing states: 100%|██████████| 1/1 [00:00<00:00, 1766.77it/s]
Linear regression... Linear regression done! (in 0.08027291297912598 sec) Running on 1 inputs (4990 steps)
sample = 500
fig = plt.figure(figsize=(15, 7))
plt.subplot(211)
plt.plot(np.arange(sample), y_pred[:sample], lw=3, label="ESN prediction")
plt.plot(np.arange(sample), y_test[:sample], linestyle="--", lw=2, label="True value")
plt.plot(np.abs(y_test[:sample] - y_pred[:sample]), label="Absolute deviation")
plt.legend()
plt.show()
Determination coefficient $R^2$ and normalized root mean square (NRMSE) :
r2_score(y_test, y_pred), nrmse(y_test, y_pred)
(0.9999869380130667, 0.0003004561614127586)
forecast = 50 # now, predict 50 steps ahead
(X_train, y_train), (X_test, y_test) = split_timeserie_for_task1(forecast)
sample = 500
fig = plt.figure(figsize=(15, 5))
plt.plot(X_train[:sample], label="Training data")
plt.plot(y_train[:sample], label="True prediction")
plt.legend()
plt.show()
sample = 500
fig = plt.figure(figsize=(15, 7))
plt.subplot(211)
plt.plot(np.arange(sample), y_pred[:sample], lw=3, label="ESN prediction")
plt.plot(np.arange(sample), y_test[:sample], linestyle="--", lw=2, label="True value")
plt.plot(np.abs(y_test[:sample] - y_pred[:sample]), label="Absolute deviation")
plt.legend()
plt.show()
Determination coefficient $R^2$ and NRMSE:
r2_score(y_test, y_pred), nrmse(y_test, y_pred)
(0.9953476933554506, 0.010759732086259792)
Let's have a look at the effect of some of the hyperparameters of the ESN.
The spectral radius is defined as the maximum eigenvalue of the reservoir matrix.
reservoir = reset_esn()
states = []
radii = [0.1, 1.25, 10.0]
for sr in radii:
W = mat_gen.generate_internal_weights(units, spectral_radius=sr,
proba=density, seed=seed)
reservoir.W = W
s = reservoir.compute_all_states([X_test[:500].reshape(-1, 1)])
states.append(s[0].T)
Computing states: 100%|██████████| 1/1 [00:00<00:00, 2545.09it/s] Computing states: 100%|██████████| 1/1 [00:00<00:00, 4198.50it/s] Computing states: 100%|██████████| 1/1 [00:00<00:00, 3045.97it/s]
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(radii)*100+10+i+1)
plt.plot(s[:, :units_nb])
plt.ylabel(f"$sr={radii[i]}$")
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()
$-$ spectral radius $\rightarrow$ stable dynamics
$+$ spectral radius $\rightarrow$ chaotic dynamics
In most cases, it should have a value around $1.0$ to ensure the echo state property (ESP): the dynamics of the reservoir should not be bound to the initial state chosen, and remains close to chaos.
This value also heavily depends on the input scaling.
The input scaling controls how the ESN interact with the inputs. It is a coefficient appliyed to the input matrix $W_{in}$.
reservoir = reset_esn()
states = []
scalings = [0.1, 1.0, 10.0]
for iss in scalings:
Win = mat_gen.generate_input_weights(units, 1, input_scaling=iss,
proba=input_connectivity, input_bias=True,
seed=seed)
reservoir.Win = Win
s = reservoir.compute_all_states([X_test[:500].reshape(-1, 1)])
states.append(s[0].T)
Computing states: 100%|██████████| 1/1 [00:00<00:00, 807.84it/s] Computing states: 100%|██████████| 1/1 [00:00<00:00, 4092.00it/s] Computing states: 100%|██████████| 1/1 [00:00<00:00, 4419.71it/s]
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(scalings)*100+10+i+1)
plt.plot(s[:, :units_nb])
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()
The input scaling can also be used to rescale the inputs and ajust their influences.
The leaking rate ($\alpha$) controls the "memory feedback" of the ESN. The ESN states are indeed computed as:
$$ s(t+1) = \underbrace{\color{red}{(1 - \alpha)} s(t)}_{\text{previous states}} + \underbrace{\color{red}\alpha f(u(t+1), s(t))}_{\text{new states}} $$where $s$ is the state, $u$ is the input data, $f$ is the ESN model function, defined as:
$$ f(u, s) = \tanh(W_{in} \cdotp u + W \cdotp s) $$$\alpha$ must be in $[0; 1]$.
reservoir = reset_esn()
states = []
rates = [0.03, 0.3, 0.99]
for lr in rates:
reservoir.lr = lr
s = reservoir.compute_all_states([X_test[:500].reshape(-1, 1)])
states.append(s[0].T)
Computing states: 100%|██████████| 1/1 [00:00<00:00, 2744.96it/s] Computing states: 100%|██████████| 1/1 [00:00<00:00, 664.71it/s] Computing states: 100%|██████████| 1/1 [00:00<00:00, 4306.27it/s]
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(rates)*100+10+i+1)
plt.plot(s[:, :units_nb])
plt.ylabel(f"$lr={rates[i]}$")
plt.xlabel(f"States ({units_nb} neurons)")
plt.show()
Let's reduce the input influence to see what is happening inside the reservoir (input scaling set to $0.1$):
units_nb = 20
plt.figure(figsize=(15, 8))
for i, s in enumerate(states):
plt.subplot(len(rates)*100+10+i+1); plt.ylabel(f"$lr={rates[i]}$")
plt.plot(s[:, :units_nb])
plt.xlabel(f"States ({units_nb} neurons)"); plt.show()
The leaking rate can be seen as the inverse of the reservoir's time contant.
During this task, the ESN is trained to make a short forecast of the timeserie (1 timestep ahead). Then, it will be asked to run on its own outputs, trying to predict its own behaviour.
units = 500
leak_rate = 0.15
spectral_radius = 0.48
input_scaling = 1.0
density = 0.1
input_connectivity = 1.0
regularization = 1e-6
seed = 1234
reservoir = reset_esn()
forecast = 1
(X_train, y_train), (X_test, y_test) = split_timeserie_for_task1(forecast)
states = reservoir.train([X_train.reshape(-1, 1)], [y_train.reshape(-1, 1)], verbose=True)
Computing states: 100%|██████████| 1/1 [00:00<00:00, 2786.91it/s]
Training on 1 inputs (20000 steps)-- wash: 0 steps
Linear regression... Linear regression done! (in 0.19173073768615723 sec)
start = 0; seed_timesteps = 100; nb_generations = 500
init_inputs = X_test[start:start+seed_timesteps]
X_t = X_test[start+seed_timesteps: start+nb_generations+seed_timesteps]
X_gen, states = reservoir.generate(nb_generations, init_inputs.reshape(-1, 1), return_init=True)
X_gen = X_gen.flatten()
plot_generation(X_gen, X_t, init_inputs, nb_generations, seed_timesteps)
Lorenz attractor is another well-known example of chaotic timeserie generator.
Here, the ESN has to predict 3 covariant values at once: the $x$, $y$ and $z$ coordinates of the Lorenz system.
We can also ask the ESN to predict one value given the two others, and so on. Of course, we can also ask for new generations of the timeserie.
from reservoirpy.datasets.chaos import lorenz
timesteps = 25000
X, t = lorenz(timesteps)
# rescale between -1 and 1
X = 2 * (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0)) - 1
plot_lorenz(X, t, 1000)
fig = plt.figure(figsize=(20, 20))
plt.subplot(211, projection="3d")
plt.plot(X_gen[:, 0], X_gen[:, 1], X_gen[:, 2], lw=3, label="ESN prediction")
plt.plot(X_t[:, 0], X_t[:, 1], X_t[:, 2], linestyle="--", lw=2, label="True value")
plt.plot([], [], ' ', label=f"$R^2 = {round(r2_score(X_t, X_gen), 4)}$")
plt.plot([], [], ' ', label=f"$NRMSE = {round(nrmse(X_t, X_gen), 4)}$")
plt.legend(); plt.show()
Not as easy as Mackey-Glass...
More fun with ReservoirPy and chaotic data: multiscroll attractor, Rabinovich-Fabrikant equations, Hénon map...
plot_attractor_examples()
import requests
import pandas as pd
data = requests.get("https://www.data.gouv.fr/fr/datasets/r/d2671c6c-c0eb-4e12-b69a-8e8f87fc224c").json()
df = pd.DataFrame(data)
df.head()
| casConfirmes | deces | date | gueris | reanimation | testsRealises | testsPositifs | hospitalises | nouvellesHospitalisations | nouvellesReanimations | decesEhpad | casConfirmesEhpad | casPossiblesEhpad | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 191 | 3 | 2020-03-02 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 1 | 212 | 4 | 2020-03-03 | 12.0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 2 | 285 | 4 | 2020-03-04 | NaN | 15.0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 3 | 423 | 7 | 2020-03-05 | NaN | 23.0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 4 | 613 | 9 | 2020-03-06 | NaN | 39.0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
French Covid-19 data is available at www.data.gouv.fr.
We selected 4 interesting features: the number of people hospitalized and new hospitalizations per day, and the number of people in intensive care (réanimation) and new admissions in intensive care per day.
plot_dataframe(feat_names, df)
plot_covid_results()